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ABSTRACT 


This thesis presents the application of EIGENVALUE 
SENSITIVITY ANALYSIS and SINGULAR VALUE ANALYSIS to the 
control of a SEA-SKIMMING supersonic missile,in the vertical 
plane. 

The study is divided in four basic parts: 

a) The development of the mcdel. 
b) EIGENVALUE SENSITIVITY ANALYSIS with respect to the 
variation of the aerodynamic parameters of the 


autopilot/airframe of the missile. 


c) Analysis of the time response with respect to the 
variation of the aerodynamic parameters. 


d) Robustness analysis and improvement of the system, 
using the SINGULAR VALUE ANALYSIS. 


All the analysis 1s based in results of simulatian 
programs using the Software available at the NAVAL 


FOSTGRADUATE SCHOOL. 
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I. INTRODUCTION 


The majority of tactical missiles used against surface 
targets are SEA-SKIMMING and fly as close as possible to the 


surface of the water in order to make difficult the 


detection and reaction of the enemy. 

This work addresses the problem of analysing the design 
of the altitude control system of a supersonic SEA-SKIMMING 
missile in both aspects of SENSITIVITY and ROBUSTNESS, with 
respect to the variations of the aerodynamics parameters. 

Those variations can develop from uncertainties in the 
model as well as from changes in the flight conditions of 
the vehicle. 

The SENSITIVITY analysis makes it possible to verify, 
for a gQGiven model which parameters are more important for 
the desired flight path. 

The ROBUSTNESS analysis has the objective of verifying 
if the system could be affected by perturbations or noise 
and will be the base on calculations of modified feedback 
Qains which are capable of producing the desired flight of 
the missile despite perturbations. 

The contents of the different sections are as following: 

Section II presents the model and the state equations of 
the system as well as the calculations to complete the chain 


Of control of the missile, considering the restrictions of 
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the SEA-SKIMMING scenario, where the altitude control system 
cannot tolerate any overshoot to avoid that the missile hits 


the water. 
Section III develops the SENSITIVITY ANALYSIS taking 


into consideration the variation of the eigenvalues and time 
response of the system of state equations for variations of 
the aerodynamic parameters in a range of +25% of the 


Original value. 

In section IV is described the ROBUSTNESS ANALYSIS based 
upon the minimum singular value of the return difference 
matrix as function of the frequency; section V introduces a 
design technique using optimization routine in order to 
obtain a "complete robust" solution and section VI presents 


the final conclusions. 
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II. MATHEMATICAL MODEL OF THE MISSILE 


The work assumes a SEA-SKIMMING missile with the flight 
in the vertical plane according to the path shown in Figure 


Ii.l. 


100.0 120.0 140.0 160.0 


ALTITUDE(FEET) 


40.0 60.0 80.0 





0.0 eo O50 We 10.0 12:5 13°59 17.5 20.0 
TIME(SEC) 


Figure II.1 - Missile Trajectory in the Vertical plane 


The missile is assumed to fly at constant velocity,at 
150 ft, when it is commanded to a final altitude of 30 ft. 


Only the pitch channel will be analysed. 
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Ae wets MISSILE 

The missile will be considered as presented by Arrow 
ERet. 11. ‘This model has been selected in order to avoid 
classification problems. The geometry as well as its 
characteristics are presented in Appendix A. 

The states that were considered in our model are defined 
in Table I as well as the relevant aerodynamics parameters. 

In the Figure II.2 we have the schematic representation 
of the uncoupled pitch control chain of the missile, 


neglecting the yaw and roll movements; 


10(s+0,. 5) 
(s+50. 0) 





Figure II.2 -— Block Diagram of the Pitch Control Chain of 
the Missile. 


where 

n - "bias" acceleration; 
zcb 

ee commanded acceleration; 
a, - achieved acceleration; 


SD 


v(t) - velocity; 
H(t) - altitude; 


H. ~- commanded altitude. 


The Transfer Function G(s) represents the missile itself 
and the other blocks are the controller. The input to this 
block will be the "commanded" acceleration and the output 
will be the achieved acceleration. 

The commanded acceleration comes from the feedback of 
the states of the controller (compensated state, velocity 
and altitude) added to a "bias" acceleration. 

In the diagram of Figure II1I.3, we have a Sepp ilies 
representation of the system with the states used in 
modelling our missile. 

The "STATE EQUATIONS" which represent our system are the 


following: 


x = -150.0 Se 20 Capes — YSIS 


xX = 1.353 Xf ee Seat + 1.334 H 
= qt Ss 


ere x, — Sx, + 6.572 9, -6.481 H, 

x= —44. 3320 | So On 

x = x, — 0.1482 x, - 0.0395 x, (II. 1) 
x, = -188.4 x_/ + 188.4 x, 
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X= -O0. 4608 i ZeLoi ea 0. 3406 ee eon ae 


“15.0949 x, - 20.13 x, - 0.143 x, + 0.4608 n,_ 
~ 0.4544 H 
Ss 
x = -50.0 x. — 495.0 x_ + 495.0 H 
3 g Ss 
cree ao 
X47 ~17-644 x. - 4.705 x, 


This system of equations can be represented in matrix 


form, where we have: 


x =AX+ Bu 

-150 O O O 2,646.64 -705.74 9O O O O 
Me Sal O 1.0 O O O O O O O 
=e al O §=5.0 O O O O O O O 
O O O O —-44.33 aye il O O O O 

= O O O i130 =§-O.410 —-0.04 O O O O 
O O O 0 O -188.4 188.4 0 O O 
eee ae ee 2 et —-20.13 -0.14 O O O 
O O O O O O O- 3065-4935 0 

O O O O O O O O O 1 

O O Oo O -17.64 Ae 7 O O O O 


17 


O O 
ses BSS 1.534 
6.572 —-6.481 


O O 
O O 
B = O O 
O O 


0.461 -0.454 


Oo 4935.0 
0 0 
u=-—F x 
0 O 0 0 O 0 O 0.049 0.98615 1.404 
a 
O O 0 0 O 0 0 0 O O 
The vector input “u" has two components 
Tee = the commanded acceleration; and 
ne - the desired altitude. 
The input 22 is considered as two parts; one 
corresponds to the feedback of the states (from the 


controller) and another is the "bias" acceleration of lg. 
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TABLE I 
DESCRIPTION OF THE SYSTEM STATES AND 
AERODYNAMIC COEFFICIENTS 


Ke eee eee ween wwe filtered acceleration 
Ke eet twee eee ween ween acceleration compensator 
Xp tte age Te a! wo sleMeiel ow ene acceleration compensator 
Ky eee ees SMenelela\a/clate ai'a ses ses Q - pitch angular rate 
Ke ee we eee eee teens 2-2 «* —- angle of attack 
Kt eee ee eee eee eee a Pitch tail incidence 
Xt eee Som commanded pitch tail incidence 
Kg eee ee eee oe controller compensator 
Ky tet eee eens es eee wan cess Sooo ge H(t) —- altitude 
Kg Tete eee eens eee ss VED £- Velocity 


AERODYNAMIC COEFFICIENTS 


ee, slope of curve of pitching moment coefficient 


on on in Cc. per degree pitch control incidence 


Cet OPE of curve of normal force coefficient Cue a 


OND change in Cy per degree pitch control incidence 


i? 


90>2°21 


Lesost— EEO 


(0 “05+S) 
3 °0+S)0T 


























(295/pes)b 
oe x 
: €rt"o 
»Ped} (224) cs°9 
do t 4 (ieee) (cist—} 
"Xx 











(1+ 2)s 


Crl‘o 
(L+ =) (cge‘o-} 

















Figure II.3 — The Missile Model 
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The model can be separated in two main parts; the 


CONTROLLER and the AUTOPILOT/AIRFRAME. 


B. THE CONTROLLER 

The controller was developed based in the procedure 
introduced by Dowdle CRef.2]. 

The missile was considered as a point mass and, 
therefore the altitude as function of time is calculated by 
double integration of the achieved acceleration. 

The controller was developed without take in 
consideration the missile itself, i.e., the transfer 
function G(s) in Figure II.2 was assumed as "unity". 

According to the CRef.2], this procedure is valid if the 
bandwidth of G(S) is larger than the crossover frequency of 
the controller. 

Due to the constraints of the scenario, the controller 
has to be designed in such way that no overshoot will be 
tolerate. Further requirements are that the angle-of-attack 
and the acceleration, have to remain below a certain level. 

The controller are responsible for the commanded 
acceleration by feedback of its states. 

The technique used to solve the feedback problem was the 
LINEAR QUADRATIC STATE FEEDBACK REGULATOR. This solution has 
a large stability margin with respect to phase variations, 


as well infinite gain margin. 
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The "cost function" to be minimized is the following: 


@m 


J = { x oily Q 2b) es ule) R u(t) + dt ieee) 
0 


In Figure I1.4 we have the schematic representation of 


the CONTROLLER. 


H 
S 
a, (t) V(t) H(t) 
1 1 + 
= _ o es 
= S Cc 


Figure II.4 -— The Controller 


The LEAD compensator indicated in the controller is used 
In order to improve the transient response characteristics. 
Its transfer function is : 


K# (Sstw_ ) 
z 


G_¢s) = 
(stw ) 
p 
where K = 10.0 
w= O.5 rad/sec 


w = 50.0 rad/sec 


In order to have the state equations of the CONTROLLER, 


we draw in FIGURE I1I.5 its signal flow graph using the 


adequate transformation of the compensator and showing the 


states that will be feedback. 


H 
s 


+ 


= 
nt) 7 i/s v(t) l/s H(t) 7 7 


eg ee eee 


a 
F ae. 
4 


cn el. Es 
Fy 1 
Z 1 x 1/s Sx wu WwW k 
i 1 Z p 
—W 
p 
1 


Figure II.5S Signal Flow Graph of the Controller 


The state vector is defined as: 


x(t) = C x(t) (H(t) -H_) v(t) J] t CT. 2) 


In the Figure II.5 the states of the controller are 


represented as indicated in Table I. 


From the diagram, we have the following STATE EQUATIONS: 


x =Ax+Bu (II.3) 
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The state matrices are: 


w Kiwoi- w_) O 
p z 
A = O O 1 
oO O O 
O 
B = O 
1 
D = C ik O J 


Since we consider to feedback all the states, let in 


equation II.1 


Q = 1 (Identity)? 


and R 


il 
© 
© 
> 


Using the OPTSYS program to solve that OPTIMUM CONTROL 


problem, the following feedback gains were obtained; 


ae = 0.049 
Fo = 0.98615 
FK = 1.4044 


correspondents to the feedback of the states indicated 


in Figure II.S. 
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C. THE AUTOPILOT/AIRFRAME 

The controller was developed assuming the transfer 
function of the Autopilot/Airframe as unity. This assumption 
implies that the achieved acceleration is equal to the 
commanded acceleration; however, they are different due to 
the aerodynamic of the missile.In order to take that in 
consideration we have to introduce the autopilot/airframe 
model. 

Normally, the sea-skimming is’) a SKID-TO-TURN missile 
(turns instantaneously), but due to problems of 
Classification we have used for our Sea--Skimming model the 
Autopilot given by Arrows CRef.1] for a BANK-TO--TURN 
missile (has to bank in order to turn) that was designed to 
fly at 30,000.0 ft instead of sea level. 

Take in consideration that the main objective of this 
work is to establish a procedure for sensitivity and 
robustness analysis and that we are studying only the 
vertical movement of the missile uncoupled of any other 
movement ,we consider that the approximation is useful. 

In Figure II.& we have the block diagram of the missile, 
representing the plant without the controller. 

The state equations derived from this diagram are on the 


Appendix A. 
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Figure II.& Block Diagram of the Airframe/Autapilot 
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D. THE TIME RESPONSE 

Using the software CONTROLS, the system represented by 
equations (II.1) was simulated and the time response for the 
altitude, velocity, acceleration and angle of attack were 
obtained. 

The altitude is in Figure II.7 and we can see that the 
missile arrives at the desired altitude in about 10 seconds. 

No overshoot is shown and a small steady-state error is 
presented characterizing our Type O system. The steady-state 
error was less than 5%. 

The angle-of-attack ,reproduced in Figure II.8, varies 
between 3.7 deg and -1.2 deg satifying the condition 
required by the airframe, as indicated in CRefs. 1 and 2]. 

The acceleration is in Figure II.9 and the velocity is 
in Figure II.10 representing a smooth behavior of the 
system. 

In the next chapter will be presented the Sensitivity 


Analysis using the developed model. 


27 


ORIGINAL SYSTEM 


meme oe te eee eeeeecemesr eee st eeu nn nen eteeeraoeeee*eeeesteoe meseneneseweeeeeeeeseeeseeeseeseeeeeeeeomeseeesercemeseeeeereensenese 


eoeoee eee Meee eseeeeBeeeeeeeeoFtseeeeeeaseeseoee ees ee Peeeee eee meses seeseHeeeeeeeeteeeaeeHeEPHet Heese Peewee weeseMesesneeeetet tone 





meme m eee ese Beweeceecen eer ren eeeFeeeuntr eee eseseees SBeeoeteenrt an MeoeeeereeeBeeeeeeeeteeeeseeeePeeeeeune ua meee srseeeseMeeeeereeB es wee 
. 


se momo ee Mme ee eee nn eeeeeee et eeeeseeeeFeeeeee ee Sunes te see Meee tees ePeeest tee SFeseteee Petes t ten Meee seeseseeseseeeeeteseesonen 


seme ec eee Puce ee te Feo eMetCeeneeeFee BF ZZ BSB SHH EHHT HSE HEHEHE HEHEHE Mee eeZr ee FePeeFeeteHHeeeeesesHHeseeees MeeeeetesemersteseerfPoeestoneosne 


eoeeeee se wmeet ere eo em eeF Fetes tPF enn sess eseeFeostenteFeeMesFeose set eHPestesestenest oes eBeeBeesese mes sePeesemeseseeteeegtooesene 


- meee eee Pu eeeeeeeaeWeseeeeeeeteeeeeeeeFeoeneeeees PEF FF meee ccc eBeseseeeetetercaseeeseaGsewmasescs Pesceesetoecvoe™Meeseteoseeteeseesen 


O'O9TO'OSTO'OFTO'OETO' OST O'OTLO'OOT 0°06 O°O8 O'O0LA O'O9 O'OS O'OFO'OE 
(LO)ACALILIV 


Figure II.7 -— Altitude(ft) vs. Time(sec) 
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Figure II.9 — Acceleratian(ft/sec’) vs. Time(sec) 
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Figure II.10 — Velocity(ft/sec) vs. Time(sec) 
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ITI. SENSITIVITY ANALYSIS 


A. INTRODUCTION 

The objective of the SENSITIVITY ANALYSIS is to 
analyse how the behavior oaf the system is affected by 
variations of the aerodynamic parameters. Those variations 
are due to uncertainty of the design and/or problems ocurred 
during the flight. 

From the developed model we can see that the important 


Gi and C. 


parameters to our system are C_s, Che? ate Ne? 


Ma 


as 


described on TABLE I. 

The variations could be interpreted as perturbations in 
our STATE EQUATION SYSTEM (A-matrix) and the analysis, 
basically, take into consideration the change in the poles 
positions due to change in the parameters. 

As we have feedback, we replace the A matrix by the 


augmented A-matrix ¢ an =A -BF ). 


The sensitivity of a matrix where the elements are 
time-invariant,can be verified in different ways. Among 
those procedures, some are enumerated below. 

First, using the method introduced by Golub and Van Loan 


CRef.3] we can analyse the matrix to verify how it is "ill- 


conditioned"; where an "ill-conditioned" matrix means a very 


sensitive system. 


oe 


In order to accomplish the analysis, we have to 


Calculate the "condition number". It is defined as: 


K(A_) = a 
A 1 | As! [| AS I | 
where 
1 | ALI! is the norm of augmented A and 
—i1 ; 
[| Al [}| 1s the norm of inverse augmented A 


The condition mumber depends upon the considered norm 
but it can be proven that the condition of the matrix will 
be the same for all norms. 


When we use the 2-norm one has 


a | 


C(A_) 

ee 
C(A_) 

= a 


K (A) = (III.1) 


a 


where 6 cae? and 6 (A) are the maximum and minimum 


Singular value of the A-augmented matrix, respectively. 

Calculating the condition number for the different 
values of the parameters it is possible to establish how the 
sensitivity of the system changes when there are changes in 
some elements of the state matrix. 

The problem of this kind of analysis is that we can not 
verify the sensitivity of the system with respect toa 
specified element of the state matrix but only the 
sensitivity of the overall system. 

In Frank (CRef.4], we have another method where the 


analysis is done based upon the STATE SENSITIVITY EQUATIONS 


that could be characterized in time domain, frequency domain 
orin terms of a performance index. 

The basic idea is to define "sensitivity functions" and 
analyse the behavior of the system by those functions. 


The analysis starts with the system defined by: 


= £ (xgtyu, a) Crriee> 


|X. 


where: 
- x represents the state vector; 
- t is time ; 
- u is the input vector ; and 
- « represents a vactor of the different parameters. 
If e& 15 changed by A« (some perturbation) and using 


Taylor’s expansion,we have: 


A o&. Cibo) 
J 


Axa CE sce) =) ax | 
Ax. oe 
j 0 


The subscript x, indicates that the partial derivative 


1s taken at the nominal parameter values. 


The partial derivative is called the trajectory vector ( 
; .th 
A >? with respect to the j— parameter and has the same 


dimension as the state vector. 
The components of that vector are defined as the 


TRAJECTORY SENSITIVITY FUNCTION. 


<p 


SG = 


at x, 
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The equation III.4 represents the partial derivative of 


en , 
the i— state with respect to the bl parameter. Assuming we 


have n states and m parameters, it will be possible to form 
an x m matrix called TRAJECTORY SENSITIVITY MATRIX (III.5). 


The matrix is also called the JACOBIAN MATRIX. 


axa Ox a 
da, a0 
pN = ; : Ciel ..2) 
axn axn 
ax aX 
1 m 


The change in the behavior of the system will be 


determined By: 


a 
4x = A (t,x) Aa =) rn. Ac. CIII.6) 
cox = J J 


ot 


j=i 
The analysis is based in how the sensitivity function 
Changes when we change the parameters. The result could be 
visualized by plotes of dX as function of time, as states of 
the expanded system. 
Finally, we consider the EIGENVALUE SENSITIVITY 
presented by Porter & Crossley CRef.S]. 
This latter procedure was choosen to be used in the 
analysis of this thesis, where the sensitivity of the sistem 


is related to the variation of the eigenvalues of the 


oo 


augmented matrix ‘ne = A - B F) with respect to the 


variation of the parameters. 
In order to further develop our analysis, the concept of 


MODAL ANALYSIS will be introduced. 


B. THE MODAL ANALYSIS 

The concept is based in generating the input vector of a 
system by linear feedback of the state vector in such way 
that the prescribed eigenvalues are associated with the 
dynamic model of the resulting close-loop system. 


As an example, let us assume a scalar system: 


~ =a x(t) + b ult); GT Tee) 


1f u(t) = Q ====> x = x, @xP (at) 


where exp (at) defines the "mode" of the system. 
If a = 0 the system is stable and if a 32 O the system is 
unstable. 
Assuming g as the feedback gain, we have: 
u(t) =gq x(t) (III.8) 
Inserting the value of u from equation III.8 in equation 


Peles 


e 


x = (a + DB gq) xt) C1Eis 9) 
and, solving for x(t): 
x(t) = x(0O) expCla + bg)t] 


Now the mode of the system is defined by the exp(atbg). 
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When we have a state vector, the system is represented 


by 
X=Ax+t+ Bu CTiso) 
Considering 
u=- F x 
and assuming all the states as observable and 


controllable, the substituting value of u in equation 


(IT1I.10), we have 


x = (A - BF) x or x =A xX ol T leet 1 ) 


The precise nature of the free motion of the continuous 
time system following any disturbance can be described in 
terms of eigenvalues and eigenvectors of the augmented plant 


matrix (A_). 
a 
Assuming A. has n distinct eigenvalues CA pA gree ADDS 


then it also has on corresponding linearly independent 


elgenvectors u_,uU_,e-.,U such that 
i 2 Nn 


Abe =X) os 2=152,... 4H (lal facut) 


Using the calculated eigenvalues and eigenvectors we can 


write the modal matrix of ae as: 
= ie ee 
1 2 n 


where each column of U is the eigenvectors of A,” 


Based in the modal matrix U, the equation (III.12) could 


be written as 
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am UF AY Cpe. tS? 


with A as a diagonal matrix of the eigenvalues. 


The matrix A is obtained from equation (III.13): 
U ae U = A (III.14) 


In addition to the eigenproperties of A the 


a’ 
corresponding properties of its transposed matrix (A? play 


an important role in the modal analysis. 


The matrix Bc has the same eigenvalues of as but a 


different sets of eigenvectors.It can be represented as 


a be — “, ay j = 1 Sayre... A 


where ay represents the eigenvectors of As: 
In CRef.5] it is shown that 


VV. Meg Oy. oe ae (III.16) 
where 85 1s the Kronecker deltas; i1.e. 


855 = i 1 f 1 = j 


and 83 = 0 if i # j 
Thus in matrix form 
V‘'U = I (III.1i7) 


where the collumns of V are the vectors Vie 


Using the modal matrix of AS and introducing a new state 


vector y(t), we have: 


x(t) =U y(t) (III.18) 
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The derivative of both sides, gives us: 
x(t) = U y(t) (III.19) 


with x(t) = Ax (t) (III.20) 


Combining equations (I1I1.18),¢(111I.19) and (III.20), 


ae 
have: 
y =U ALU y(t) 
yY a oA 
and using the equation (III.14), 
y(t) =A y(t) fret) 
Because the matrix A is diagonal the equations are 


uncoupled and the solutions of the system represented by 


(III.16) are given by 


y, ©) = y, €0) exp (rat) 


(III.22) 
Aplying equation I11.18, 
x(t) = Cu, Us see uid Vis 
Ye Clete 25) 
Yn 
which implies that 
= + eee 
x (t) uUsy, oo exp (a,t) + Un¥5 10) exp (rat) 


Taking the value of y(t) in equation (IT1I.18), 


yt) = Ux ct) (LIT. 24) 
Using equation (III.17), 
y,; 6) = V5 x (oO) CIIT.25) 
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Combining equations (III.23) and (III.25): 


x (t) - x (O) uUsVvs exp (r,t) CLL e226) 


Equation III.26 shows that the system represented by 


(III.8) thas its modes described by the eigenvectors of A 
and its transpose and eigenvalues of AS matrix. Also, in 


order to have the system stable, x(t) ---> 0 as t ---> @®@ , 


which happens if and only if 


Re ds < O , for all the values of i 


In the Modal Control CRef.3Sj all this procedure is shown 
for the case when some of the eigenvalues are equal which 
does not apply to our system. 

As the role of the eigenvalues is so important, we can 
verify how the sensitivity of our controll system by 
analysing the change in the eigenvalues of the matrix A or 
the augmented A for the close loop system when any element 
is changed. 

In order ta avoid calculating the new eigenvalues for 
each different A, the equation (III.26) can be used to find 


approximated values that can be analysed analytically. 
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Ge EIGENVALUE SENSITIVITY 

In our system we consider the change in the elements of 
the matrix due to variations in the aerodynamic 
coefficients. 

Basically, we concerned with the eigenvalues that are 
changed in the direction of the right hand plane,making the 
system unstable or affecting the time response. 

The result of the analysis could be plots or tables 
relating the eigenvalues with the variation of the 
parameters. 

For each change in the parameters, we have a new matrix. 

Approximations to the eigenvalues associated with the 
different plant matrices in the vicinity of the original one 
could be used without having to compute the eigenvalues of 
each new matrix. 

We describe the approach presented in the reference 
CRef.5] applying the "first order eigenvalue sensitivity” to 
the change of one of the parameters and doing the complete 
analysis with the calculated eigenvalues. 


1. First-order Eigenvalue Sensitivity 


From equation III.11, we have: 
x(t) = A_ x(t) 
a 


let’s assume 


*o 7 | a] 
a kl 
where aud is any element of As: 
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Differentiating equation (IITI.12) with respect to 





any indicates that 
aA uu. + A sui = @dAl u. + A dul CL Lila 7, 
aa x da aa = 3a 
kl kl kl kl 


multiplying (III.27) by Me then gives 


vi @Au. + vi A aul =v: @dAl uu. + A. Vv. Bul 
‘i da - ga aa : . ; aa 
kl KI kl KI 





As we Ais equal to di Vi » we have 


Vv: @AU, = Vv: BAI U 
a2 1 


Fi kl 





Which reduces to the set of equations 


SAi =v. (k) uj, (1) 5 isk=l=1,2,...,n (III.28) 


J4.4 


Due to 


aa 


day 


ll 
a 
Gy 


In equation (III.28), Vv; tk) and u, (1) represent the 


you element of V3 and the ycls element of Us respectively. 


Those coefficients may be considered as the 
elements of a set of n eigenvalue sensitivity matrices. 


The new eigenvalues (a5) of the altered matrix, 


when the element and is changed will be calculated as 
A = xr + : ri Aa 
1 1 Sau Kk] 
kl 
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When we have more than one parameter altered, the 
linearity could be used and the final eigenvalue will 
obtained from the summation of the sensitivity coefficients 


multiplied by the respective change in the coefficients. 


As an example of this procedure, let’s consider our 


system, changing the value of ko that 


A = 0.02 


Bes 


When Tap is changed as specified, the following 


elements of the augmented-A matrix of the state equations 
presented in the Chapter II are changed: 


A,(4,6): Aa 14.778 
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A_(7,6): Aa 5.035 


7 
The eigenvalues and eigenvectors of the augmented A 
matrix and of its tranpose were calculated using the program 


Matlab. 
a) Eigenvalue A, | 


oes Sloe 70 + Jeo oo 


A A, SEN A) u,(6) 4] A Au, + Cv (7) u,(6)4 A ass 


Vanier) = -0.0148 - j 0.0039 


ll 
— 


v,67) 


u, (6) = 0.014 - j 9.9027 


dA, = -159.815 + j 18.85 
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b) Eigenvalue A, | 


Te —-159.815 - 3 18.85 (complex conjugate of A,) 


c) Eigenvalue A, | 


A, = -50.0 

OA, = Cv (4) u,(6)] 4a, + Cv,(7) u,(6)] aa_, 
v_(4) = 0.0003 

v7) = BOL OGs7 

u,(6) = -0.001 

A. = -49.99998 


d) Eigenvalue Ay 


AY = -8.6i + 3 7.91 

Ad, = yt) u, (6) ] 4a, + Cv (7) u, (6) ] cverre 
Vv, (ae = -0.16467 - j 0.1322 

Vv, 67) = i 

u, (6) = -0.0334 - j3 0.132 

x = -8.95 + j 7.6355 


iy 


e) Eigenvalue A, F 
r = -8.95 — j] @eesoa 


s 


#) E1genvalue A, 3 


44 


A = -1.82 + j3 2.02 

Ad, = Cv, 4) u, (6) ] ar, + Cv, (7) u, (6)] has, 
te =e es ey Ono rot 

Viet) = -0.5737 + 3 0.1689 


u, (6) = 0.01 + 3 0.0363 
» = -1.846 + j 2.05 
gq) Eigenvalue A, * 


ne = -1.846 - j 2.05 


h) Eigenvalue A, 5 
d — ae Ge | 
AA = Ev (4) u_(6)] Aa + Ev (7) u_(&)] Aa 
8 8 46 8 8 7 


v,(*) = 0.1741 


YZ) = -0. 6597 
u, (6) = —-0.0396 
d = -2.681 


d = -0.9 


3 


Cv (4) u (6) J Bay * Cv (7) uO) ] Aas, 


< 
+ 
ll 


0.1718 + j 0.016 
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v,(7) = -0.5737 - j 0.1689 
u,(6) = 0.01 - j 0.0363 
d = -0.4741 


9 


j) Eigenvalue Dion 


D,o = ENA k a? uo) Sar, + Eve, 67) u,, 64] 4a, 


Vv, (4)= 0.1741 
10 


Vv. (7)= -0.6597 
a) 


U,, (d= -0.0396 


dae = -—0.1104 


For the analysed parameter , Our system i5 
unsensitive, taking in consideration the small differences 
between the new and original eigenvalues. 

Repeating this procedure for all different parameter 
values will give us the approximation of the new eigenvalues 
and that can be helpfull in a qualitative analysis of what 
Darameter has more effect over our model. 

AS we want to verify the effect in the time 
response and due to the Software available the analysis will 


be done completely with the calculated exact values. 
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2 Exact Eigenvalues 


Using the Controls Program, we will verify the variation 
of the eigenvalues of the augmented A matrix ( A — FB) with 
respect to the variation of the parameters and by simulation 
we analyse the effects on the time response. 

At the end of the analysis we compare the results of the 
example used for the first order approximation with the 
exact values that were calculated. 


As reference, the eigenvalues of the original system 


ares 
A, = -159.746 + 3 18.979 
h, = -159.746 - j 18.979 
B= 47992 
A, = -8.60697 + j 7.90968 
A. = -8.60697 - § 7.90968 
RK, = -1.81854 + j 2.0238 
A, = -1.81854 - j 2.0238 
AR = =2. 70266 
A, = -0.503819 
A. = -0.143073 
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Each aerodynamic coefficient, that takes part of the 


considered model was changed up to +25% and the eigenvalues 


of the perturbed system matrix (A-augmented) were analysed 
as well as the TIME RESPONSE of the system . 


Ae C 
ma 


On TABLE II we have the eigenvalues related to the 


considered values of Cd as we can see, with exception 


of the Fe este yale the changes are small and the system 
is pratically insensitive to that variation. 

The time responses are shown in Figures III.1, III.2,I1II.3 
and III.4 . 


- Angle-of-attack (Figure III.1) —- changes in the same 
direction of Cw with the minimum value changing about 


O.5 degrees in the negative direction. The maximum 
Value stays almost constant ; 


- Acceleration (Figure III.2) -—- the change was less than 


2 
2 ft/sec either the maximum or the minimum values; 
- Velocity (Figure III.3) - almost the same 3; and 


- Altitude (Figure III.4) - almost the same. 


b. C 
m§p 


The change in the eigenvalues are on TABLE III 3 the 
time response 1s on Figure III.S,II1I.6,I1III.7 and I1I1I.8. 
The variation of the eigenvalues are slightly higher 


with changes in Si Dae the system stays stable and the time 


response has a small variation. 
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Angle-of-Attack (Figure III.S) —- The minimum values 
increases negativelly about 0.5 degrees but in the 
contrary direction of the coefficient, i-e., increasing 
the value of Snes negativelly, the minimum was smaller; 


Acceleration (Figure III.6) - small changes; 
Velocity (Figure III.7) - small changes; 


Altitude (Figure III.8) - pratically the same. 


ewok 


No 


The change in the eigenvalues are on TABLE IV and the 


time response is on Figures III.9,III.10,III.11 and III.12. 


on 


The changes in the eigenvalues are small but the effect 


the time response is more sensible, mainly on the angle- 


of-attack. 


time 


Angle-of-Attack (Figure III.9) - the maximum and the 
minimum increase 1.0 degree in the positive = and 
negative directions, respectivelly; 


Acceleration (Figure III.10) —- maximum changes of 102 
in the negative direction (minimum) and 20% in the 
positive direction (maximum).These values are smaller 


when we increase the value of Ce : 


Velocity (Figure III.11) - small changes; 


Altitude (Figure III.12) - does not change. 


d. 7e 


N Sp 


The change in the eigenvalues are on Table V and the 


response are on the Figures I11.13,I1I1II.14,111.15 and 


Dia. 


Angle-of-Attack (Figure III.13) - Changes almost 1.0 
degree with the maximum and minimum increasing when the 
Value of the parameter increases; 
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Acceleration (Figure III.14) - the maximum increases 
with the value of the parameter ; 


Velocity (Figure I[I1I.15) - pratically the same for all 
the considered values of the parameter; 


Altitude (Figure III.1i6) - constant. 


eee 


TABLE II 


C - SENSITIVITY 


=). OS 





-159.741 + j 18.9843 -159.7435 + j 18.9814 
-159.741 - j 18.9843 -159.743 - j 18.9816 
-49.9992 -49.9992 | 


someaee2 tog 6.357130 —8.563586 + j 7.18736 


poesdgess ~=" 967371350 —-8.56386 - 3 7.18736 


—-2.31740 + j3 2.95574 —=2.14117 + j3°2.432/75 





—2.31740 —- j 2.95574 Seeneii7 seams 4s275 

-1.79944 —2.14200 

-0.351650 —-0.510282 | 

-O. 14304 -90.14306 | 
a, -0.07 -0.08 

seove746 +9 1629762 eave, ol ee 1eao7S5 

eis7.74959 -ej he. 9762 SAD7ef al “ej Terr 7 so | 

gee 7 7 7 2 =89 7 9992 | 


—3.66197 + j 8.559389 -8.72178 + 3 9.14021 


-8.66197 + j 8.955389 -gevai./o + eee. 14021 


{ 
-1.47981 + j 1.81088 -1.22077 + j 1.70081 | 
-1.47981 - j 1.81088 -1.22077. - j 1.70081 
S020 1 -3.67125 
| -0.49779 -0.49148 
| -0.14308 -—0.14308 


a 


o1 


ate 


-162. 604 


-162. 604 


-49.9992 


—3. 746035 


—3. 74603 


-1.24694 


-1.24694 


—-3. 86370 


-0.49103 


-0.14311 


m §p 


—-158. 243 


=138. 243 


-49.9991 


-10.1621 


-10.1621 


-1.47981 


-1.47981 


—2.2/254 


-0.50786 


-0. 14307 


TABLE III 


=O05 07 


ad 


=Oore 


Casp7 SENSITIVITY 

-0.06 
+ j 16.7131 -161.200 
= ten 1S 1 -161.200 
-49.9992 
+ j 8.86120 -7.13033 
- j 8.86120 =7VSOSS 
+ j 1.66681 -1.54230 
- j 1.66681 -1.54230 
-3.30548 
-0. 49828 
-0. 14309 

-0.09 
+ j 19.9192 -156.678 
- 5 19.9192 -156.478 
-49.9990 
+ j 6.98960 ~11.7920 
- j 6.98960 -11.7920 
+ j 1.981088 -2.03752 
- j 1.981088 -2.03752 
-2.02276 
-0.51120 
-0.14306 


+ 


+ 


+ 


fe eS 1 


17.9141 


8.50992 
GE S0992Z 
1.803335 


12 S80sSs 


20.7492 


20.7492 


S/aoal7 
ae. 2eo 19 
2.00411 


2.0411 


RR St ER Rn SS SSS SS SSS ei IG 


TABLE IV 


C —- SENSITIVITY 








-159.750 + j 18.9895 -159.748 + j 18.9842 
-159.750 - j 18.9895 -159.748 - j 18.9842 
-49.9992 -49.9992 
-8.76667 + j 8.27012 ~8.68784 + j 8.09494 
-8.76667 - j 8.27012 -8.68784 - j 8.09494 
-1.18629 + j 1.74321 -~1.47163 + j 1.84903 
~1.18629 - j 1.74321 -1.47163 - j 1.84903 
-3. 63456 S227 75 

-0. 48878 -0. 49663 

-0.14314 -0.14310 

men. 0.16 O.17 

| 159.744 + j 18.9736 -159.741 + j 18.9683 
! -159.744 - j 18.9736 -159.741 - j 18.9683 
| -49.9992 -49.9992 
-§.52441 + j 7.71274 -8.44077. + j 7.50213 
—8.52441 - j 7.71274 -8.44077. - j 7.50213 
=-2.15317 + j 2.32249 -2.40105 + j 2.67918 
(220153517 ~ “@ees2249 -2.40105 - j 2.67918 
| -2.20610 -1.88565 

-0.51041 -0.51651 

~0.14305 -0.14303 | 


23 





CN sp 
-159.004 
-159.004 
-49.9996 
-9.37545 
-9.37545 
-2.24432 
-2. 24432 
-1.98723 
-0.51356 

ee 


CN sp 


-160.106 
=al6é0S106 
—49.9990 


-8. 25204 


Se. re ef gr 
- 


—-8. 25204 


aii57570 


-1.57370 


-3.18783 


-0. 49848 


-0.143190 


12.1138 


i222 So 


ToS? 


Pa 19 


2.00109 


2.26109 


21.9632 


21.5632 


8.25703 


8.25703 


1.85345 


1.85345 


TABLE V 
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C = SENSETIViITY 
sp 


0.03 


=159: 748 


—-159. 748 


-—49.9992 


-8. 68784 


—-8.68784 


—-1.47163 


-1.47163 


—3 22277 5 


—-0. 49663 


-0.14310 


Sad 


+ 


+ 


0.06 


-160. 459 


—-160. 459 


-49.9988 


-7.91449 


-7.91449 


155009 


—-1.35009 


—~—3.60970 


-0.49275 


-0.14312 


+ 


+ 


+ 


18.9842 


18.9842 


) 8.09494 


8.09494 
1.84903 


1.84903 


235.6487 


25.8437 


8.57438 
8.57438 
1.75196 


1.75196 


D. CONCLUSION 

AS a conclusion for the sensitivity analysis, we have 
the following: 

With respect to the eigenvalues, the pitch control 


system 15S more sensitive to changes in Cc. and CN Sp than 


changes i1n Ce and Cues ; and 


Considering the time response, we can verify that the 


angle-of-attack 1S sensitive to all the parameters with 


higher Variation with Cie? the acceleration 1S more 
sensitive to C and C ; the velocity has the maximum 
Max Na 


value changed less than 10%2% with respect to all =the 
parameters; and the altitude is insensible for the 
considered variations. 

It is not necessary to emphasize the importance of this 
analysis for the design of the control system, but thinking 
only in the eigenvalues and comparing the results of the 


first order approximation presented for a iee with those of 


fapre IfF for C changed to -0.06 (ac = 0.02) we can have 
m&p m &p 


an idea of the importance of the parameter, using that 


approximation. 
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IV. ROBUSTNESS ANALYSIS 


A. INTRODUCTION 

In the previous chapter we take into consideration the 
change in the eigenvalues and in the time response of the 
system due to variations on the aerodynamic parameters that 
take part into the model; under the robustness analysis we 
verify what we can expect of the system when it is affected 
by any kind of perturbation. 

The main idea is ta use the robustness analysis to 
design a new set of feedback gains in order to have the same 
behavior in different situations and/or environments which 
act as perturbations to our missile. 

In classical frequency domain, for single-input single- 
Output systems(SISQ), a robust design can be achieved using 
Bode, Nyquist or Nichols plots. With these techniques it is 
possible to define gains that give us gain and phase margins 
for a "robust" system. For multi-input multi-output systems 
(MIMO), the classical techniques are no longer valid. 

Taking this in consideration, we have analysed the 
robustness of the original design by means of the minimum 
Singular value of the return difference matrix for different 
frequencies and we will use the same technique to improve 


the control system. 


72 


In order to show why we are applying this technique, we 
start with the Nyquist Criterion for a SISO system and 
extend the analysis to a MIMO system, using the minimum 
singular value theory, as explained by Lehtomaki, Sandell 


and Athans in CRef.6]. 


B. THE NYQUIST CRITERION 


Let's take in consideration the SISO system represented 


imerrguse I[V.1. 


ee) | Gis) can ee 


Figure IV.1 Single-Input Single-Output System 


where G(s) includes the controller and plant dynamics 
with unity feedback . 

The Nyquist criterion states that if the open loop 
transfer function G(s) does not have any pole in the right 
half s plane, then the locus of G(s) will not encircle the 
point (-1,9) in the Nyquist plot, where jw is substituted 
for s and the axes represent the real (Re) and imaginary 


(Im) parts of Gljw) for various frequencies. 
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ae 


A Nyquist Plot 1s given in Figure IV.2 and it can be 
used to design feedback gains that ensure a robust system. 
The measure of the system stability are Gain and Phase 


margin obtained from the Figure, as indicated. 


__— _ perturbed 


GM system 


Figure IV.2 Nyquist plot —- Stable System 


The Phase Margin (PM) is the additional phase lag at the 
gain crossover frequency required to bring the system to the 
limit of instability and the Gain Margin (GM is the 


reciprocal of the magnitude [Gljw) | at the frequency where 
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the phase angle is -180°. The gain crossover frequency is 
the frequency at which |G(jw) | is unity. 

The Phase and Gain margins are a measure of how close a 
polar plot is to the -i + jO point. 

Considering the system subject to additive perturbation, 


as we have in the diagram of Figure IV.3. 


Figure IV.3 Additive Perturbation 


In order to have the system stable, despite the 
perturbation, we need 
lAGCjwd) | < |1 + GCjw) | (IV.1) 
As we can see in Figure IV.2, this condition ensures a 
stable system. 
This idea can be extended to MIMO problems through the 
use of matrix norms and of applying the Multivariable 


Nyquist Theorem. 
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C. MULTI-INPUT MULTI-OUTPUT SYSTEM 

The generalization of the SISO theory discussed in the 
previous section has been made for the MIMO problem. 

Let’s consider our system as represented in chapter III 
and apply the transformation in order to have the open loop 


transfer function between input and output. 


x =Ax + Bu GIVE) 
Applying Laplace Transform to the equation (IV.2), we 
have 


sX = AX + BU (IV. 3) 
where xX and U represent the Laplace Transform of qx and 


u(t), respectively. Taking the value of X in equation (IV.1) 


1 


X = CsI - A] BU (IV. 4) 
As Y=cC X (IV.5) 
We have Y =C {Csi - Al 'B (IV. 4) 
U 


This transfer function corresponds’ to G(s) in Figure 


IV.4. 
u t—| em St 





Figure IV.4 Return Difference 
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From Figure IV.4 it is possible to define the "Return 
Difference Matrix" (Kwakernaak and Sivan, (CRef.6]) that 
corresponds to opening the loop at a ~<given point and 
connecting an external input variable at this points; for 
example, at point 1. Assuming the input u as zero, we have: 


Y,¢s) = - G(s) H(s) V (5s) (IV.5) 
where V (s) corresponds to the new input variable and 
>? the "returned variable". 
The difference Detween Y,¢=s) and V(s) is 
Vis) - Y,¢s) = € I + GlsdH(s)] Vis) (IV. 6) 


The matrix I + G6(s)H(s) 1s defined as the "Returned 
Difference Matrix", with G(s) equal to 


Be ee eee (1V.7) 


As explained by Lehtomaki ,Sandell and Athans in CRef.7], 
the Multivariable Nyquist Theorem is’ derived from the 


relationship 


detCI + G(sd)H(S)] = Gel (s) CIN 9) 
€o1 (s) 
where $c 1 (s) = det(sI - A) corresponds to the 


characteristic polynomial of the open loop transfer function 
and Fol (s) = det(sI - A + BF) corresponds to the 
Characteristic polynomial of the closed loop’ transfer 


function from the system represented in Figure IV.4. 
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The multivariable Nyquist theorem which requires that a 
closed loop stable system have the same number of 
counterclockwise encirclements of the origin by the locus of 
the determinant of I + G(s)H(s) as the number of open loop 
poles that are unstable. 

If I + G(sdIH(sS) 165 quasi-singular, a small change in G 
may make the matrix singular ; this causes the detCI + 
Gt(sd)H(s)] to become zero and the Nyquist encirclement count 
to change indicating an unstable system. 

Basically,in order to analyse the robustness of the 
system we have to verify how close the return difference 
matrix is to being singular as a function of frequency (w). 

The natural measurement of the singularity 1s the 
minimum “singular value", since this is the tightest norm. 

The singular value of a matrix M 1s defined as 


6 = ALC mim) 327? (IV. 10) 


mt represents the conjugate transpose matrix of a 
generic M matrix and A; any elgenvalue of the product of M 


times M. 

If the minimum singular value is close to zero the 
matrix 1S qQuasi-singular and the system is not rabust. 

If one assumes 


6(I + GljwoHCjw)) > Ky (IV.11) 
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then it can be shown CRef. 7] that the gain and phase 


margins for the system may be represented by equations 


@WWS12) amd (I1V.13), i. e@., 


Gain margins: GM = (IV.12) 


+ cos Ci - «0'] (IV.13) 


= 


it 


Phase margin: PM 


Universal gain and phase margins curves, presented by 
Gordon CRef. 81] and repeated in Figure I[IV.5 allow the 
designer to pick a singular value that corresponds to a 
desired stability margin for the system to be designed. 

From singular value plots the we will be able to 
identify the critical frequencies or range of frequencies 
where the singular value is below the chosen level from the 


Figure IV.93. 
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Figure IV.5 - Universal Gain and Phase Singular Value Plot 
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D. SINGULAR VALUE ANALYSIS 

Using the program developed by Gordon CRef.8], with 
small changes to make it compatible with the size of our 
system and to input the matrices more easily, we calculate 
the minimum singular value of the return difference matrix, 
with the frequency varying from 90 to SO rad/sec. 

The program utilizes the IMSL subroutines to generate 
the return difference matrix as well as to made the singular 
value decomposition, given the minimum singular value for 
different frequencies. 

The return difference matrix depends on the point where 
we open the loop represented in Figure IV.5. For the 
explained case, considering point 1, we Rave the "output" 
difference matrix, that corresponds to I + G(s)H(s)3 if the 
loop 1s open at point 2 we have the “input” case and the 
return difference matrix will be I + H(s)G(s). 

From the Figure IV.5 we choose as reference, for our 


analysis, the minimum singular value as 0.6 that gives us a 


phase margin of about 35° to characterize a "robust" system. 
The minimum input singular value is plotted in Figure 
VO 
Lt shows that the system is robust, except for 


frequencies between 0.6 and 5.0 rad/sec. 
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The results related to the minimum output singular value 
are in Figure IV.7 and the robustness of the system is poor 
from O to 18 rad/sec. 

Those results indicate that the behavior of the system 
under perturbations should present problem for frequencies 
below 10 rad/sec. 

To confirm this expectation, we have plotted the BODE 
diagram (magnitude) for the closed loop system considering 
the angle-of-attack and altitude with respect to both inputs 
(acceleration and desired altitude). 

The diagrams are on Figures IV.8 to IV.11, where can see 
the attenuation at the low frequencies indicating that we 
can not expect good response at those frequencies and one 
could further expect problems with robustness as is evident 
from the singular value plots. 

On the next chapter we try to improve the system 
designing a new set of feedback gains, in order to have a 


robust system. 
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V. IMPROVING THE DESIGN 


Under this chapter we try to improve the system by 
calculating new feedback gains that will yield a more robust 


system, taking into consideration the results indicated in 


Figures IV.6 and IV.7. 

The POPLAR program, developed by Gordon (CRef.8] will be 
used. A numerical optimization technique is applied in order 
to increase the minimum singular value of the considered 
return difference matrix, therefore, a resultant robustness 
of the design. 

The development in this chapter will be preceded by 
comments on numerical optimization and a description of the 


computer program. 


A. OPTIMIZATION 

The optimization was accomplished by means of the 
Automated Design Synthesis Program (ADS) developed by 
Vanderplaats (CRef.9]. 

The purpose of ADS as of others numerical optimization 
routines is toa find the "best" possible solution for the 
problem, starting from an initial set of variables and 
updating the design iteractively. The problems can develop 


in convergence of the method and in the computer time 


89 


needed. he the problem has multiple solutions, the 
optimization does not always lead to the absolute optimum. 

The ADS program is’ designed as a Black box optimizer 
which aiwees the user to choose combinations of one 
dimensional search, optimization algorithm and strategy. 

ADS is used as a subroutine and the parameters that 
correspond to the different applications are chosen Dy the 
user, as explained in CRef.91]. 

The preferred method used in the ADS applications is 
referred to as Sequential Unconstrained Minimization 
Techniques and can be considered as a method that starts 
with an objective function and the constrains combined into 
an augmented objective function and then minimizing this 
function as if the problem was unconstrained. 

ADS employs penalty function techniques as well as an 
Augmented Lagrange Multiplier. 

As presented in CRef.8], the iteration between the 
user ’s program and the ADS routine can be represented by the 
block diagram of Figure V.13; where the "info" parameter is 
used as a "flag" for a dialogue between the user and the ADS 
routines. 

The details for utilization of ADS are presented by 


Vanderplats in CRef.8]. 
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Figure V.1 Organization of ADS Program 
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B. POPLAR PROGRAM 

The Gordon's program has the capability of calculating 
the feedback gains that improve the response of the system 
either By pole placement and/or by increasing the minimum 
Singular value of the return difference matrix. 

The main objective is to improve the robustness by 
increasing the minimum singular value above a desired level 
chosen from the universal curve reproduced in Figure 1V.5. 

"Optimum" values of the feedback are calculated by 


minimization of an Objective function whose pole placement 


part is 
1 
OBJ =) ( ag - Ag + CAL = AY )* 
di i di 1 
1=1 
where 
Ae - real part of the desired eigenvalue i; 
di 
he - real part of the computed eigenvalue; 
1 
a imaginary part of the desired eigenvalue 1; 
di 
Ay imaginary part of the computed eigenvalue i; 


and, for the minimum singular value optimization, we 
a 2 
OBJ 2. £ max CO, 6, — & (jw,p) J +; 
J 


where oy indicates the desired minimum singular value 


and 6 1s the minimum singular value at a certain frequency. 
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The optimization procedure changes the feedback gains 
until the minimum singular value is raised to the desired 
level. 

The pole placement and robustness Program calculates the 
return difference matrix to be considered, working in the 
complex space (as the analysis is conducted in the frequency 
domain) and the objective functions, calling the ADS routine 
to make the optimization. 

Input’s for the program are the matrices of the state 
representation of our missile ( A,B,C and F),where the 
observation matrix (C) was specified as Identity because we 
are assuming all the states are observable as well as 
controllable. 

Other inputs are the desired minimum singular value, 
desired pole locations and frequency interval to be 
considered. 

In our application of POPLAR the initial values of the 
feedback gains are those of the original system and the 
parameters for application of the ADS program are: 

Strategy —- Augmented Lagrange multiplier; 


Optimizer - BFGS variable metric method for unconstrained 
minimization; 


One-Dimensional Search —- Polynomial interpolation. 
Two situations will be considered; the INPUT —- minimum 
Singular value, where the return difference matrix is I + HG 
and QUTPUT-minimum singular value with I + GH as the return 


difference matrix. 


9S 


C. INPUT -— MINIMUM SINGULAR VALUE 

The program was first run to improve the minimum 
Singular value for the input situation within the same 
interval of frequencies used in Chapter IV. 

Due to the large CPU time involved, the program was run 
for the critical Band of frequencies and the result 
evaluated for the complete interval. 


The minimum singular value taken as reference from the 


Universal curve was 0.6 that implies a phase margin of 55% 

We have an objective function taking into consideration 
the improvement of the minimum singular value as well as a 
part corresponding to the pole placement. The main goal is 
to improve the singular value but with a solution such that 
the poles are in positions where the time response satisfies 
the requirements. 

The weight of the pole placement part was considered as 
10% of the singular value part. 

The program was started considering only part of the 
feedback matrix as free parameters, increasing this number 
until a reasonable solution was found. 

The best solution was found keeping the original 
feedback gains and calculating gains to feedback all the 
other states to generate the commanded acceleration. 

The computer output of the best solution with the 


corresponding inputs is given in Appendix B. 
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The feedback gains that improve the robustness of the 


system are the following 


: | S2O% —1. lameness —le7-7707 O2Ge 7.07°0-049 0.986 1.4 
; Oo O 0 O 0 Oo O Oo O 0 | 
The plot of the minimum singular value for different 
frequencies iS given in the Figure V.2, along with the 
values for the original system. 
AS we can see the minimum singular value was increased 
in the critical values of frequency with the minimum 


changing from 0.40 to 0.63 that indicates an improvement in 


the phase margin from 20° to 35°. Also the bandwidth where 
the system was less robust was changed from 5.0 rad/sec to 
approximately 1.0 rad/sec. 

The robustness of the control system for all the 
considered frequencies is better, arriving close to 0.9 with 
a considerable improvement at low frequencies. 

Using the software Controls, the time response was 
plotted in Figure V.3. 

The maximum and minimum values of the angle-of-attack 
are smaller than those of the original system; the maximum 
acceleration and velocity are slightly lower and =€_4n 
undershoot of about 10% appears on the altitude, that could 


be a problem in some flight conditions. 
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We now turn to a consideration of the output return 
difference matrix although this is less important from a 


practical viewpoint. 
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Figure V.3 - Time Response for the Improved System (INPUT) 
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D. OQUTPUT — MINIMUM SINGULAR VALUE 

In the Figure V.4 we have plotted the output minimum 
singular values that correspond to the run where the input 
singular values were improved, as well as the original 
values. 

As a considerable degradation of the output singular 
values was evident, we try, in the last part of this work, 
Ze) arrive to a situation where the system could be 
reasonable robust in both cases, increasing the output 
singular value even with some reduction at the input 
Situation. 

After several runs of the program, using different 
combinations of optimization techniques with different 
starting points, ist was verified that due to the 
characteristics of the plant, that increases in the input 
singular values resulted in decreases in the output 
Situation. 

In order to obtain some improvement, an effort was made 
to modify the system. 

The first step was to look at the controllability matrix 
of our model, represented by 

1.855632D-01 4.948168D+02 
2.725664D+13 —-2.690449D+13 
2.725664D+13 —-2.690449D+13 
2.861121D+00 —-2.822009D+00 


2.201811D+01 2.171710D+01 
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—6&. 78177 7p 2eo 6.689111D+00 
-2.286975D+00 2. 255563D+00 
9.855433D+00 -9.719614D+00 
—~1.508382D+01 1.487685D+01 
4.905487D-01 -4.264398D-Ol1 

As we can verify, all the elements are non ~ero, this 
means that the poles can placed at any desired position; but 
a close look at the 2nd and 3th rows indicates very high 
order numbers compared to the others. This”) situation 
generates some numerical difficulties in pole placement. 

Using engineering judgement, we have to define what kind 
of changes have to be made for achieving a more robust 
system. 

First we try to get a balanced A-matrix by changing 
units, i. e., the angles will be in radians instead of 
degrees and the pitch angular rate in radians per second 
instead of degrees per second. Despite some reduction on the 
numbers was not possible tO arrive to a satisfactory 
solution. 

More positive effect was obtained by adding the effect 
of the angle-of-attack in the controller (see Figure II.3). 

As the original system has the output singular values 
higher than O.5 for frequencies above 10 rad/sec, we have 
kept the original feedback gains and calculated gains to 
feedback oanly the angle-of-attack, pitch rate and commanded 


actuator. 
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The gains are the following: 


O O O 0.037 0.027 060 -0.02 0.049 0.986 1.404 


O 0 0 O O O O O O O 


The result is presented in Figure V.5, compared with the 
Original singular values and those from the improvement in 
the input case. 

The output singular values, compared with the original 
was Slightly improved for frequencies below 5.0 rad/sec. 

Resultant input singular values are plotted in Figure 
V.6, and they are lower than those from the improvement in 
the input situation but a little higher than the original. 

The time response were practically unchanged with 
respect to the original model and it is presented in Figure 
ely 

Further analysis is needed to obtain a solution with 


high output singular values. 
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Figure V.46& Final Result in the Input Singular Value 
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Figure V.7 Final Time Response 
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VI. CONCLUSIONS 


The sensitivity analysis based in the change in 
eigenvalues represents an important role in the design. 

As was demonstrated for the studied missile its behavior 
corresponds to a system insensitive to the variation of the 
considered aerodynamic parameters. 

The analysis gives a good indication of what parameters 
have to be precisely determined. 

The robustness analysis also demonstrated how useful the 
singular value analysis 1S aS an auxiliary tool for the 
designer. 

Singular value analysis indicates that the system is 
robust to input perturbations but is deficient in robustness 
to output perturbations. 

The physical nature of the problem indicates that the 
prime concern should be with respect to the input 
perturbations; the output perturbations are of minor concern 
but should be kept in mind if unusual conditions should be 
encountered by the missile. 

On the MIMO design a commomly used method is the Linear 
Quadratic analysis where the performance levels are reached 
by adjusting weighting terms in the "cost" function, but the 
results of that method, for non-diagonal R matrices, do not 


necessarily imply a robust system. Using the return 
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difference matrix we can improve the robustness of the 
design. 

One advantage of this procedure is that it permits a 
high level of interaction between the system and the 
designer. 

Further problems of a computer nature (i.e., large CPU 
time) are encountered with a high number of states. Analysis 
must be confined to a small frequency range of low singular 
Values and the complete system return matrix singular value 
are calculated with the determined feedback gains. 

One aspect that should be considered in future 
development is to improve the cpu time used when handle a 


system with a high number of states. 
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APPENDIX A 
MISSILE DATA 
Under this Appendix the missile sizing, mass properties 
and aerodynamic parameters are presented as given by Arrow 
C[Ref.1], including the figures. 
The missile is 1/6 scale of the actual circular missile 
configuration and 18S reproduced in Figure A.1. It is tail 


controlled using four identical control surfaces located 


o 
with +30 dihedral. 


A. GEOMETRY AND MASS PROPERTIES 

In table VI we have the size and mass properties, with 
the respactive values, used for development of the state 
equations. 

Only the uncoupled pitch channel was considered, 


assuming no roll movement. 





Table VI 


MISSILE GEOMETRY AND MASS PROPERTIES 


Weight w 29525 lbs 
Length 1 168 in 
Diameter d 24 in 
Reference Area for Coeff. = T et 
2 
Moment of Inertia about y I 804 slug-ft 
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Maximum body 
diameter, 4.0 


Moment reference location 


All dimensions 
in inches 


Figure A.1 Circular Configuration of the Missile 
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OT 


B. AUTOPILOT/AIRFRAME 

Figure A.2 shows the block diagram of the autopilot used 
to develop the model of the pitch control channel of the 
missile. 


The state equations that correspond to the autopilot are 


the following: 


2) = oO scat soso ne 
x = 1.353 x +x. - 1.353 9 
2 i 3 Le 
x= -6.5728 — 59.0 co. + 6. 5me ae 
z 1 3 2 C 
x = =44.3535196 *  — 592511 x 
4 5 6 
x =o 0.1482 x « + 0.0395 a 
Rea -188.4 os + 188.4 x5 
x = -0.4608 «9 —eeees! x = 0 Seog est 
7 i 2 3 4 


="lase72 we 20.13 ‘a O. 1430 cos 0.4608 ee 


The input for the autopilot is the commanded 


acceleration --- 90. 
ZG 


The values related to the aerodynamic are the following: 


Ce = - 0.06 
m sp = - 0.08 
CN = O.15 
Nf = 0.04 


Poasanee re Ssure - GQ == ———-————— 1650 lb/ft’ 


Velocity ee 3825.46 ft/sec 


Tespoy DtweudAg 





(29s/pes)b 


40 ReNQ2AY 


(pea) 
ad 
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Crt o 


=e 





ite 


(b+ 29) (ors1-) 











(r+ 2)s 
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Block Diagram of the Autopilot/Airframe 


Figure A.2 
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APPENDIX B 
POPLAR PROGRAM 
The computer output where the improvement in the input 


Singular value was found is listed below. 


PiLEs LATSPROB OUT Al 


PAGE 00) 
1 
<<<<¢¢ RUN B® $0901 >>>>> 
WMAX s @.70000 INIT FREQs @.70000 FREQ STEPS 6.10060 
WEIOHT 1s @.1C00OMELOHT 20 O.SO0COMEIGHT Js @.90000 
SYMINI« 6.60006 SVHINOs 0.80000 Rds 6.10000 1DGs 9 
ses THE A PLANT MATRIX e808 . 


-156.00000 6.060000 @.00000 9.00009 suanessacae795 78003 0.00000 6.009000 
6.900008 6.006000 
1.35316 6.00000 1.00000 @.90008 6.00000 6.00000 6.90009 8.069000 
0.90009 6.00609 
-6.57200 0.96000 -$.00000 0.900000 6.000006 @.006098 6.90000 6.96009 
0.90008 @.00000 
9.96090 9.00009 9.00000 6.00000 -44.33162 -S9.11000 6.00000 6.00006 
9.00009 0.00000 
9.00000 9.00009 0.06000 1.00000 <@.14020 -=0.039590 @.00000 6.90000 
6.00006 0.000090 
0.99606 9.000080 0.00060 0.900006 O.0000G-100. 59999 100. 39999 @.00000 
0.00009 0.00060 
0.46006 =2.25000 <8.34058 2.23000 —15.09850 -20.1539009 -0.14300 0.00000 
6.90008 @.eccce 
6.60600 @.90060 6.60006 6.90008 6.60008 6.00008 @.00008 -50.800090 
-495.00008 0.96000 
9.00080 9.00000 0.90000 6.00000 @,.00000 6.00000 6.909008 6.00000 
6.60009 1.00006 
6.90000 @.90000 6.00068 @.00000 -17.64400 —4.70508 6.09008 6.00000 
6.90006 9.90009 


pee THE B CONTROL INPUT MATRIX co 


-90008 @.00008 
333519 1.33446 
57200 <6.40100 
09000 @.00009 
00000 @.90009 
600090 6.60008 
46000 0.454596 
o0008 495.90000 
eoce8 0.00000 
oocee 0.90008 


wean THE C OBSERVATION MATRIX #88 


1.00009 6.00609 6.906000 6.00008 @.00000 @.00000 6.06008 6.00000 
@.6000090 6.00008 
6.99000 1.00060 @.00006 6.906008 @.00008 @.00000 6.00000 @.90000 
9.00000 6.00600 
0.00009 @.00008 1.90000 9.06000 6.00000 6.00000 6.00008 @.90000 
6.00609 @.00000 
6.00000 @.00000 9.00008 1.606008 @.06000 @.09000 6.00600 0.60008 
6.00090 @.000¢60 
@.eeece @,.0e00080 @.90008 @.ceeee 1.00000 @.00008 @.09000 @,.e8006 
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PILE: LATSFROE OUT 


6.00008 
0.00000 
0.00000 
6.00008 
9.00000 
8.00000 
0.00008 
6.00000 
1.00000 
8.00008 
8.00008 


ann THE 


0.00000 
0.98615 
0.00000 
0.00000 


@.80000 
6.0000¢ 
0.00000 
6.00000 
6.00008 
@.00000 
6.00000 
0.00000 
8.00000 
@.00000 
1.00000 


Al 


6.00000 


@.000ce 


0.00800 


@.e00dce 


FP PEEDOACK MATRIX ene 


@.00000 
1.40440 
@.00000 
9.00000 


9.00000 


@.c00000 


0.00000 


@.00008 


@.00000 


0.00008 


ane THE ORGEREG COMPLEX EIGENVALUES CINPUT) 


-159.80008 
-~159.60000 
=-30,00000 
=-8.00000 
=-8.00000 
=7.00000 
=-$.00000 
-1.00000 
-1.00000 
-t.00000 


=19.89999 
19.69999 
0.00000 
7.90008 
7.90008 
0.00000 
@.00000 
@.0000e 
6.00000 
6.80006 


ane THE INITIAL DATA 15 ene 


nee THE ORDERED COMPUTED EIGENVALUES 


-1$9.74232 
=139.74232 
-30.02896 
-8.42668 
6.61668 
2.46024 
-1.86206 
-1.86206 
=0.30480 
“0.34313 


-16.97508 
18.97508 
@.00000 
-7.8902)3 
7.89021 
@.co0ce 
=-2.158078 
2.15078 
0.00000 
@.00000 


eee THE ORDERED ELGENVECTORG 


-29.46394 
-4.43356 
3.60236 
8.26904 
-1.0¢6860 
@.13223 
=-1.30795 
$.29377 
=68.9027) 
-0.04065 
7.7372? 
0.09812 
8.00034 


-14,65166 -29.46394 14.65180 
~21.39685 8.13620 0.90006 
9.00008 3.79569 @.eccce 
@.09535 0.26904 0.09535 
1.56396 2.19106 @.e000e 
9.90000 0.00672 0.90006 
-0.4634) 1.30793 0.46341 
=6,64100 —11.29992 @.00008 
@.0000@ -01.51242 @.00080 
=0.56439 0.04065 0.36439 
1.42416 1.37471 0.00008 
@.00000 -6.00145 6.000868 
6.00197 @.00054 -@.86197 


1.cooce 


6.00006 


@.90000 


@.00000 


0.00000 


1.36392 
4.47925 


O.32167 
0.04739 


-1,.$6693 
-4,42362 


9.28462 
O.72477 


@.00551 


1.00008 


0.00000 


6.00000 


@.ccc00 


@.e0000 


@.e0cce 
-1.32631 


2.900090 
-4,25205 


0.99000 
0.53061 


@.00008 


9.04900 


=4.43350 
4.47925 


-1.0¢6880 
0.04739 


3.29377 
=4@.42302 


@7,73727 
@.72477 


@.36706 
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11.39655 
1.32631 


~1.36596 
9.08130 


4.64100 
4.23203 


~1.40416 
6.53061 


=8.36376 


PACE G02 


FILE: LATSPROS OUT Al PAGE 003 


0.86700 0.36376 -0.S88603 9.00000 -0.31322 $.08299 0.31522 =-0,08299 
0.28464 0.00000 -0.26068 @.00000 
-0.01491 -1.00106 -0.01491 1.00106 -0.24803 @.00000 -1.36366 =0.96010 
-1.36566 0.96010 0.49679 0.00000 0.23844 -0.01917 0.25044 0.01917 
0.19194 @.00000 0.20211 0.00000 
-O.10317 -0.18008 -0.10317 0.18008 -O.17996 0.00000 -1.34360 0.86683 
-1.34360 0.86653 0.49030 0.00000 0.25586 -0.02171 9.25506 8.02171 
0.19143 0.00000 0.20199 @.00000 
0.06046 -0.00966 -0.06046 0.00966 <251.18878 0.00000 =1.20696 =1.435904 
~1.28696 1.43904 -12.89149 9.00000 2.72214 0.84036 2.72214 =-§$ .84836 
~-140.76167 0.00000 -1828.26782 0.00000 
0.00004 0.00018 0.00004 -0.00018 0.00042 0.00000 @.07780 0.03327 
0.07750 -0.03327 1.32212 0.90000 =-0.24893 -0.$1431 0.248593 0.$1431 
14.08348 0.00000 184.98924 0.00000 
“0.00381 -0.02884 -0.00501 0.02004 -0.02110 0.90000 =0.40679 =-0.09804 
-0.40679 0.89684 =3.25274 9.00000 -0.64822 1.40662 0.64822 1.40662 
-7.11066 9.00000 =26.47836 0.00000 
eae THE F FEEDGACK MATRIX e458 


0.00000 0.00000 0.00000 0.00000 9.00000 0.900000 @.90000 0.04900 
0.98615 1.40440 
0.00000 0.00000 0.01000 0.00000 0.00000 @.00000 0.09000 0.00000 
0.00000 0.90000 


FREDUENTY MIN ADD IN SV 2ND ADD IN SV MAX ADD IN SV MIN PROD ADD IN SV 
0.70000 0.88284 3.28329 3.28329 1.341106 
FRED SvaDMo SVADXD SVHI1M SVM1X SVHNO SVMXD 


0.70000 0.05214 35.16096 0.$22202221.19727 2.00000 2.00000 


AAABA pooooD $$S38SS 
A a oO o.6 CUS 

A A oO db. Os 
AAAAAAA D D $3538 
A A D D $ 
A A D 5 $ 
A A pooooo $3$33s 


FORTRAN PRODRAM 


FOR 


AUTOMATED DESION SYNTHESIS 


VERSION 1.00 


CONTROL PARAMETERS 


PILGs LATSPAOD OvT Al 


ISYRAT @ 
1GAAD e 


OBJECTIVE FUNCTION VALUE 





OPTIMIZATION RESULTS 





DESIGN VARIABLES 


@.10670F °06 


1PRINT = 10006 


LOWER UPPER 
VARIABLE BOUND VALUE BOUND 
1 -0.10000E°0S5 O.3S9SI7E-O1° 0. 1000008 
2 =0.10000E*90$ <0. 11494890) 0.30000E°05 
3 -0.10000EF +05 0.30785E-0) @.10000E°0S 
4 =0.10000F°0S -O.169SS3E°01 @.10000E°0S5 
5 ~0.10000E°OS -0.90698E°O1 @.10000F°0S5 
é =-0.10000E°05 0.31634E°00 @.10000F°0S 
? -0.10000E 05 8. 7074SE°0) @.10000E°08 
PUNCTION EVALUATIONS « 233 
eee THE OPTIMIZER COMPUTED OUTPUT 15 sen 
see THE ORCERED COMPUTED EIGENVALUES 
-1S9.73711 -19.06007 
-159.73731 19.04807 
-$0.03259 0.00000 
6.04686 -7.97487 
-6 04668 7.97467 
-7.06468 0.00000 
-5.,00031 8.00000 
0.37434 -0.39758 
-0.37434 0.39758 
0.14305 0.00000 
eee THE ORDERED EIGENVECTORS 
~27.97833 -12.15036 -27.97833 12.13036 1.71545 9.00006 
~0.61352 0.96862 2.93255 @.00000 2.22260 @.6000006 
6.46230 0.57794 -3.69572 0.00009 
0.27374 0.08283 0.27374 0.08283 0.34835 @.00000 
0.168090 9.03811 14.59894 6.00009 -2.89061 6.00009 
0.03296 0.00429 0.00719 @.00008 
“1.33070 0.40259 -1.335070 0.40259 -1.69691 6.00000 
-O.91973 0.16036 <+-72.37257 6.00006 14.44669 6.00008 
-0.1897% 0.01996 01.48660 @.90000 
-0.06026 -8,34608 -6.06026 6.34808 -6.30621 6.90000 
0.50392 @.49214 1.91362 @.00008 -6.60346 6.00006 
0.02459 -6.00252 0.00666 0.00000 
0.00060 0.00187 0.00060 -0.00187 0.00597 6.00000 
=0.00022 0.06230 0.27422 6.00006 @.17743 0.00008 
-0.03261 -0.04087 8.26187 @.00000 
-0.04631 -9.96230 -0.04631 0.96239 0.26536 0.00006 
9.13517 0.04574 0.43442 6.00000 -0.20900 6.00000 
0.02460 0.03047 -0.19580 0.00000 
-0,10883 6.14151 0.10053 @.14151 -6.19409 6.00006 
0.12746 @.04950 0.41611 6.00000 -6.20337 0.90006 
0.02446 0.01046 -0.19600 @.00000 
70,06095 -6.00637 -0.06093 @.00637 -243.063660 6.000086 
0.19754 6.04027 @.F9126 6.00008 -0.67280 6.00008 
20.2269 -14.34287 1762.69556 6.00008 
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-0.61352 
6.462350 


6.16800 
0.03296 


-0.91973 
-0.10979 


6.30392 
@.02459 


=0.00022 
0.03261 


0.13517 
6.02460 


@.12746 
0.02440 


O.19754 
20.20692 


-0.96862 
-0.57794 


-0.03511 
6.00429 


0.19036 
O.01796 


-6.49214 
6.60252 


6.06230 
6.94007 


-0.04576 
-0.03047 


=6,64956 
-0.93046 


=-0.04027 
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